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ABSTRACT 

We examine the combined effects of winds and photoionizing radiation from 0-type 
stars on embedded stellar clusters formed in model turbulent molecular clouds 
covering a range of masses and radii. We find that feedback is able to increase the 
quantities of dense gas present, but decreases the rate and efficiency of the conversion 
of gas to stars relative to control simulations in which feedback is absent. Star 
formation in these calculations often proceeds at a rate substantially slower than 
the freefall rate in the dense gas. This decoupling is due to the weakening of, and 
expulsion of gas from, the deepest parts of the clouds’ potential wells where most of 
the star formation occurs in the control simulations. This results in large fractions 
of the stellar populations in the feedback simulation becoming dissociated from 
dense gas. However, where star formation does occur in both control and feedback 
simulations, it does so in dense gas, so the correlation between star formation activity 
and dense gas is preserved. 

The overall dynamical effects of feedback on the clusters are minimal, with 
only small fraction of stars becoming unbound, despite large quantities of gas being 
expelled from some clouds. This owes to the settling of the stars into virialised and 
stellar-dominated configurations before the onset of feedback. By contrast, the effects 
of feedback on the observable properties of the clusters - their U-, B- and V-band 
magnitudes - are strong and sudden. The timescales on which the clusters become 
visible and unobscured are short compared with the timescales which the clouds are 
actually destroyed. 
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1 INTRODUCTION 

Star formation occurs inside giant molecular clouds (GMCs) 
in an hierarchical fashion. At the largest scales, clouds 
convert a small fraction, typically a few percent, of their 
gas reservoir to stars before dispersing on timescales of 
~ lOMyr. Most stellar groupings dissolve into the field 
on a similar timescale, but several classes of identifiable 
objects, such as open clusters, OB associations, scaled OB 
associations and globular clusters, survive for much longer 
times. What it is that governs which of these paths a given 
young stellar population takes is not clear. 

There are many processes acting on different timescales 
which can disperse a gro up of stars into t he field, recently 
reviewed and discussed bv iKruiisserJ (|2012| ). The mechanism 
most intimately connected with the star formation process 
itself, and acting on the shortest timescales, is referred to 
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as ‘infant mortality’ or ‘infant weight-loss’ (iLada fc Ladal 
l2003l : iBastian fc GoodwinI l2006l) . These authors gathered 
together observations of embedded clusters and showed that 
the rate at which embedded clusters are formed is 10-20 
times higher than what would be expected from the rate at 
which gas-free open clusters are formed, if all embedded 
clusters were to evolve into open clusters. They therefore 
inferred that ~90% of embedded clusters do not evolve 
to become open clusters, but are instead dispersed on the 
same timescale (5-lOMyr) as that on which they evolve 
from the embedded to the exposed or open phase. This 
implies that the expulsion of gas that makes the clusters 

visi ble is also what d estroy s them. _ _ 

iFall et ahl (l2005l ) and IWhitmore et ahl (l2007l ) use the 
large number of clusters visible in the Antennae galaxies to 
plot mass-age diagrams in which they also deduced that 
~ 90% of clusters are lost in each logarithmic age bin, but 
that the cluster mass function was unaffected. The dispersal 
mechanisms acting on each timescale seem therefore to be 
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mass-independent. _ 

This issue was revisited bv iGieles fc BastianI (120081 ') who 
instead examined the most massive clusters in logarithmic 
age bins in the Antennae, M51, the SMC and LMC, M33 
and M83. Th eir results for the Ant ennae were in agreemen t 
with those of iFall et al.l (120051 ') and lWhitmore et al.l ll2007f) . 
and they also found evidence of mass-independent cluster 
disruption in M51, but not in the other systems studied. 
iBastian et al.l ll2009ll corrected for the non-constancy of 
star formation in the Antennae merger system and found 
evidence for cluster disruption taking place on timescales 
< lOMyr, but not on longer timescales. The observational 
picture of clusters at very young ages is thus somewhat 
murky. 

From the theoretical point of view, there are two 
potential ways of ensuring that a recently-formed cluster 
disperses on a similar timescale to that on which star 
formation goes to completion and the cluster becomes 
exposed. The first is gas expulsion and this has been 
extensively investigated. Most such studies have relied on 
analytic models or N-body simulations where the gas is 
represented, in one way or another, as an ex ternal potential 
which is removed over some t i mescale (e.g Tutukoy 1978; 
Hills! [1985: iLada et al.l 1 19841: iB oilv fc Kroupa 2003a h; 


Goodwin fc BastianI ^Od: Baumgardt &: Kroupa 12007 : 

Pfalzner fc KaczmarekI 2013lh These simulations have 


generally found that gas removal while the potential is 
still gas-dominated unbinds large fractions of the stars. 
However, hydrodynamic simulations call into question 
whether the potential i n which the stars are situated is 
ever gas-dominated (e.g. Offner et al.ll200^ : iKruiissen et al.l 
l2012l: iGirichidis et al.ll2012ll . 

With improvements in algorithms and increases in 
computing power, it has become possible to model both 
star formation and stellar feedback self-consistently in hy¬ 
drodynamic simulations. Some authors have concent rated 
on th e effects of therm al accretion feedback (e.g iBatel 
l2009bl : lUrban et al.l 2010|) or jet s /outfiows from low -mass 
stars (e .g. |Li fc Nakamura!|2006l: Cumim^hnm^t_aT IMl 


Krumholz et aki ' 20121 : Federrath et ah! 20141 : Mvers et ah! 

2014^ . These works have generally concentrated on the 


effects of feedback on setting the masses of individual stars, 
the determination of the stellar initial mass function, and 
the setting of the overall star formation efficiency. These 
feedback processes are expected to have limited influence 
the large-scale dynamics of the clouds or the clusters, 
except insofar as they reduce the rate at which gas is 
converted to stars. 

At larger scales, modelling of the destructive effects 
of HII regions on whole GMCs has been popular, since an 
HIT region can in principle expand and clear out a large 
fraction of the volume of a molecular cloud on a relative 
short timescale. These simulations have produced somewhat 
mixed result s, depe nding on the clouds under investigation. 
!Walch et al.l (l2012f) modelled the effect of ionising radia¬ 
tion from central O-stars on clouds with various fractal 
dimensions and find it to be highly destructive, dispersing 
their 6.4pc radius, 10‘^Mq clouds in a few Myr. Similar 
disruption timesc ales in clou d s form ed by coll iding flows 
were reported by ! Colin et ah! (l2013! l. However. !Dale et ah! 
ll2005li . who modelled ionising radiation emanating from an 
O-star at the hub of a network of filamentary accretion 


flows, !Peters et al.l (l2O10! l (modelling irradiation from 
O-stars forming near the centre of a massive disk) and 
!Dale fc Bonnell! (1201 l! l (simulating the influence of massive 
stars on a 10®Mq cloud with a high escape velocity) all 
observed their HII regions to flicker as they were swamped 
by neutral gas delivered by accretion flows. In these calcu¬ 
lations, the dynamical influence of photoionisat ion on the 
scales of whole clouds was m uch more modest. !Dale et ^ 
(120121) and !Dale et ah! (l2013bl) examined this issue in a suite 
of simulations spanning a GMC mass-radius parameter 
space ranging from IO'^-IO^Mq in mass and 5-180 pc in 
radius. They confirmed that accretion flows restrict the 
ability of HII regions to disrupt clouds. However, they also 
showed that the clouds’ es cape velocities were crucially 
important (as predicted by !MatzneJ (l2002l) for example), 
since HII regions cannot expand at speeds much exceeding 
the (roughly constant) sound speed in gas photoinised by 
O-s tars of « 10km s~^. _ 

!Walch fc NaabI (l2014h investigated the combined 
influence of photoionisa tion and supernov a explosions, 
building on the work of !Walch et ah! (l2012l) . They found 
that the effect of the photoionisation phase before the 
supernova explosion was to delay somewhat the transition 
of the supernova remnant from the Sedov-Taylor phase 
to the radiative phase, by reducing the density of the gas 
encountered by the remnant. This allows the supernova to 
deposit « 50% more momentum into the cold gas. 

Instead _ of _ photoionisation, 

[Pelupessv fc Portegies Zwart! (|2012|) model the effects 
of winds and supernova explosions from massive stars on 
an embedded cluster. They do not model star formation, 
but they do investigate the influence of feedback on an 
admixture of stars and gas. They find that the influence of 
feedback depends very strongly on the efficiency with which 
the clouds retain the injected energy (which they choose to 
parameterise), but that star formation efficiencies as low as 
5% can result in bound systems surviving gas expulsion. 

The second possible explanation for the apparent poor 
survival ch a nces of embedded clusters, investigated by 
!Clark et ah! (12005! ). is that the clusters are never bound in 
the first place. This seems plausible on the fa ce of it, since 
GMC s exist with a variety of virial ratios fe.g !Dobbs et al.l 
[2^. However, two problems exist. Firstly, star formation 
in unbou nd clouds tends to result in flat stellar mass 
functions (IClark et al.!!200^ ). Secondly, since the stars are 
able to decouple from the gas dynamics at early ages, the 
cloud being unbound does not guarantee that the stars will 
be so. 

As a counterpoint to the foregoing discussion, there has 
been a resurgence of the idea that the space distribution 
of stars should more properly be t hought of as hierachical, 
rather than merely clus t ered (e.g Elmegreen fc ElmegreenI 


!2001| : !Bonnell et ai.![200^ : !Bastian et al.!!2007l ). This is partly 

driven by the difficu lty of defining observ ationally what is 
and is not a cluster. !Bressert et ah! (|2010!) highlighted this 
particularly strongly in their survey of the surface density 
of nearby YSO’s in which they were unable to identify any 
distinc tive scales wh i ch co uld be used for such a definition. 
While !Gieles et al.l (l2012l) showed that this result does 
not necessarily imply that stars are not formed in bound 
clusters, it does show that observa tionally defini n g and 
identifying clusters is non-trivial. !Bastian et al.l (l2007l) 
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obtained similar results at much larger size-scales in their 
study of M33, in which they were unable to identify any 
preferred size scale for young stellar systems. We do not 
discuss our simulations in this context (Parker & Dale, 
2015, in preparation). 

This paper forms part of a series of studies of the 
effects of photoionisation and/or winds from 0-type stars 
on a parameter space of model GMCs constructed to reflect 
th e gross properties of the Milky Way clouds documented 
in iHever et al.l ll2009ll . In the hrst six papers, we performed 
controlled experiments involving ionisation or winds with 
the obj ect of disen t anglin g their individual effects. 

In IPale et ^ ll20l4l . we combined the two forms of 
feedback and computed quantities global to our model 
clouds, such as the star formation efficiency, average star 
formation rate and unbound gas mass, to see how they 
varied across the parameter space. The results of the study 
can be briefly summarised as follows: The additional effect 
of winds on the dynamics of the cold gas was minimal, 
although the structure of the ionised gas was strongly 
altered in many of the simulations, being compressed by 
the winds into a thin shell lining the inner walls of the 
feedback-blown bubbles. The overall influence of feedback 
was strongly dependent on the clouds’ escape velocities, 
with the lower-mass lO^M© clouds having large fractions of 
their gas reserves unbound or expelled, whereas the 10®M© 
clouds were much less severely damaged. Star formation 
rates and efficiencies followed a similar trend, being reduced 
by factors up to ~ 2 in the low-mass clouds, but scarcely 
changing in the larger objects. Feedback left the clouds 
substantially permeable to ionised gas, ionising photons 
and supernova debris. 

We previously focussed largely on the gas content of 
the clouds, mentioning the stars only as sources of feedback 
and sinks of gas. In this paper, we redress the balance 
by considering in more detail the properties of the stellar 
populations and clusters formed by the clouds. We examine 
the interplay of stars and gas in more detail and compare 
our results with more nuanced observational diagnostics 
than blunt instruments such as global star formation 
efficiencies. 

This paper proceeds as follows: Section 2 contains a 
brief recapitulation of our numerical methods. Section 3 
describes the results derived from our simulations, and our 
discussion and conclusions follow in Sections 4 and 5. 


2 NUMERICAL METHODS 

We have investigated the formation of embedded clus¬ 
ters under the influence of stellar feedback using a set of 
Smoothed Particle Hydrodynamics (SPH) simulations of 
turbulent GMCs covering a mass-radius parameter space 
running from 10‘*-10®Mq in mass and 2.5-180pc in radius. 
The clouds initially have shallow Gaussian density profiles, 
with the density contrast between the centres and edges of 
the clouds being approximately three. The clouds were given 
divergence-free turbulent velocity fields satisfying the scal¬ 
ing relation P(k) oc k~‘^, scaled to give them initial virial 
ratios of either 0.7 (which we refer to as ‘bound clouds’) or 
2.3 (which we refer to as ‘unbound clouds’). 


The therm al properties o f the neutral gas are governed 
by a piecewise iLarsonl ll2005l ') equation of state defined by 
P = kp^ with a density-dependent adiabatic exponent given 
by 


7 =0.75; 

7 = 1.0; pi < 

7 = 1.4; p 2 < 

7 = 1 . 0 ; 


P < Pi 
P ^ P2 
P < P3 
P > P3, 


( 1 ) 


and pi = 5.5 x 10“^®g cm“®,p 2 = 5.5 x 10“^®gcm“®, pa = 
2 X 10“^®g cm“®. The initial properties of the clouds are 
given in Table [T] 

The clouds are allowed to evolve, with the turbulence 
freely decaying, until they have formed a few massive stars 
or a few subclusters massive enough to host such stars. 
Stars and subclusters are represented by sink particles. The 
mass resolution of the 10"^ M© and 3 x 10“* M© clouds are 
1 and 3 M© respectively, and sinks are taken to represent 
stars. Their accretion radii are 5 x 10~®pc and their mini¬ 
mum creation densities are set to 7 x 10^cm“®. In the 10® 
M©, 3 x 10® M© and 10® M© clouds, the mass resolutions 
are 10, 30, and lOOM©, and the sinks represent subclus¬ 
ters. Their accretion radii are set to either 0.1 or 0.25pc, 
whichever is less than 1 percent of the host cloud’s initial 
radius. The minimum creation density of these objects is set 
to 4 X 10®cm“® We then enable photoio nisation and wind 
feedb ack using the algorithms pr esented in lPale et al.l (l2007l . 
l2012lf and I Dale fc Bonnelll (120081 '). In the 10“* M© and 3 x lO"’^ 
M© clouds, sinks above 20M© are assigned an ionising lu¬ 
minosity according to 


log(gH) = 48.1 -k 0.02(M. - 20M©). 


( 2 ) 


For the more massive clouds, sinks are treated as small clus¬ 
ters with Salpeter mass functions between 0.5 and lOOM©. 
The total mass contained in stars more massive than 30M©, 
M 30 is computed and the ionising luminosity of the sinks is 
set to (M 30 / 3 O) X Q 30 s“^ with Q 30 being the ionising lumi¬ 
nosity of a 30M© star from the above equation. 

Wind momentum fluxes are set in a similar manner from 
mass loss rates given by 


M(M*) 



* 

[0.3 exp (^—j 

- 0.3 


X 10-®M©yr-\ 


(3) 


and terminal velocities by 

VooiM:,) = [10®(M* - 18)°'^'^ -k 600] km s"\ (4) 

The clouds are then evolved for a further 3Myr (or as close 
to this as possible) to establish the combined effects of 
ionisation and winds in the interval between the formation 
of the first O-stars, and the first supernova explosions. 

The glo bal effects of pr e-supernova feedback were 
discussed in iDale et ^ (l2014l l. Here, we focus more on 
the detailed properties of the stellar clusters and the 
smaller-scale interplay between stars and gas with reference 
to recent observations. 
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Run 

Mass (Mq) 

Ro(pc) 

{n{H 2 )) (cm ®) 

VRMS,o(km s 1) 

(ToKK) 

(Mo) 

tff.o (Myr) 

A 

10® 

180 

2.9 

5.0 

143 

6.5 

19.6 

B 

10® 

95 

16 

6.9 

92 

11.2 

7.50 

X 

10® 

45 

149 

9.6 

52 

20.7 

2.44 

D 

10® 

45 

15 

3.0 

92 

4.9 

7.70 

E 

10® 

21 

147 

4.6 

52 

4.9 

2.46 

F 

10® 

10 

1439 

6.7 

30 

19.0 

0.81 

I 

lO'* 

10 

136 

2.1 

52 

4.6 

2.56 

J 

lO'* 

5 

1135 

3.0 

30 

8.5 

0.90 

UZ 

10® 

45 

149 

18.2 

52 

39.3 

2.9 

UB 

3 X 10® 

45 

45 

10.0 

68 

18.9 

6.0 

UC 

3 X 10® 

21 

443 

14.6 

36 

37.9 

1.9 

UV 

10® 

21 

148 

12.2 

52 

26.3 

3.3 

UU 

10® 

10 

1371 

8.4 

26 

25.6 

1.1 

UF 

3 X 10^ 

10 

410 

6.7 

28 

19.7 

2.0 

UP 

10-* 

2.5 

9096 

7.6 

18 

27.9 

0.4 

UQ 

lO'^ 

5.0 

1137 

5.4 

30 

14.9 

1.2 


Table 1. Initial properties of clouds listed in descending order by mass. Columns are the run name, cloud mass, initial radius, initial 
RMS turbulent velocity, the initial mean gas temperature, the initial mean turbulent Mach number, and the initial cloud freefall time. 



Time(Myr) Time (Myr) Time (Myr) 


(a) Run I (b) Run J (c) Run UF 




(d) Run UP (e) Run UQ 


Figure 1. Comparison of the dense gas mass (left axis and solid lines) and star formation rates (right axis and dashed lines) in the 
control (blue) and dual—feedback (black) Runs I, J, UF, UP and UQ simulations. 
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(a) Run I (b) Run J (c) Run UF 




(d) Run UP (e) Run UQ 

Figure 2. Dense gas mass depletion times in units of the freefall time at pcrit (~ 0.35 Myr) in the control (blue) and dual—feedback 
(black) Runs I, J, UF, UP and UQ simulations. The red dashed line corresponds to a depletion time of tff(pcrit)! which would indicate 
star formation proceeding at the freefall rate in gas of this density. 



(a) Run I (b) Run J (c) Run UF (d) Run UP (e) Run UQ 

Figure 3. Barcharts recording the final fate of all gas whose density comes to exceed pcrit Runs I, J, UF, UP and UQ simulations at 
the ends of each calculation. Lefthand columns are from simulations, righthand columns from dual-feedback calculations. Black: involved 
in star formation. Red: ionised. Blue: acquires higher density. Green: acquires lower density. 


3 RESULTS 

3.1 Dense gas depletion times 

iLada et al.l l|201Ch infer a linear relationship between the 
SFR and the mass of gas above a density threshold pcrit of 
10^ cm~®. We examine this relationship in our simulations 
in Figure [T] where we plot M(pcrit) as a function of time 
as solid lines using the left-hand axis scale. M(pcrit) in the 
control simulations (blue solid lines) varies by moderate 
factors of a few over the simulation durations. In the 
feedback runs, there is generally more dense gas than in 


the companion control simulation by factors of a few, and 
again with variations over time by factors of a few. 

Star formation rates (dashed lines and right hand 
axis scale in Figure [T]) are generally higher in the control 
simulations by factors of up to three compared to the corre¬ 
sponding feedback calculation. In both sets of calculations, 
the star formation rates vary non-monotonically with time 
by factors of approximately two. 

Since stars form from dense gas, M(pcrit)/SFR can be 
thought of as the depletion time of the dense gas due to star 
formation. If there are no forces present to resist gravity 
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and gas exhaustion has not yet occurred, the depletion 
time would be expected to be close to the freefall time at 
Pcrit, tff(pcrit) « 0.35 Myr. There are no magnetic fields 
present in our calculations and only one of our simulations, 
Run F, comes anywhere near gas exhaustion, so we expect 
th e control simulatio ns to form stars at the freefall rate. 
In IPale et ^ (|2014| ). we showed that star formation in 
the control simulations indeed proceeds at close to the 
freefall rate defined by the freefall times of the entire 
clouds, whereas star formation in the feedback simulations 
was generally slower by factors up to two. In both sets of 
simulations, the star formation rates also vary with time, 
although also only by factors of approximately two over the 
timescales shown here. 

In Figure [21 we show the depletion times of gas with 
p > Pcrit (defined as the instantaneous quantity of material 
satisfying this condition divided by the instantaneous star 
formation rate) in units of tfr(pcrit) as functions of time 
in the control (blue lines), and dual-feedback (black lines) 
Runs I, J, UF, UP and UQ. The red dashed line indicates 
a depletion time of tfi(pcrit). The depletion times in the 
control simulations are within a factor of two of this value 
for the virtually the whole duration of all simulations. 

The depletion times in the feedback runs in contrast 
are substantially longer, since the star formation rates 
in these models are slower and the dense gas masses are 
generally larger. The depletion times are almost always 
at least a factor of two longer than those in the control 
runs at the same epoch. Depletion times in the feedback 
runs also vary much more strongly with time, sometimes 
being more than an order of magnitude higher than in the 
corresponding control calculation, although the variation is 
often strongly non-monotonic. 

Overall, the roughly constant depletion times in 
the control runs indicate that star formation in these 
simulations proceeds quite steadily at the freefall rate in 
the dense gas. Since tff(pcrit) is a constant, and the star 
formation rates also do not vary greatly, the invariance 
of the depletion times in the control simulations implies 
that the masses of gas above pcrit, M[pcrit) are also nearly 
constant. 

In general, the depletion times in the feedback sim¬ 
ulations are longer than tfi(pcrit), sometimes by factors 
approaching ten, indicating a decoupling of the star forma¬ 
tion rate from the freefall time in the dense gas. There are 
two extremal possible explanations for this; star formation 
at densities at or above pcrit (i) proceeds efficiently but more 
slowly (ii) proceeds at the same rate but less efficiently. 
In reality, the explanation may be a combination of these 
factors. 


3.2 Production of dense gas 

We investigate the issue of the long dense gas depletion 
times in the feedback simulations by first examining how 
much dense gas the simulations actually produce. We 
identify, in all the above simulations, every SPH particle 
whose density ever exceeds pcrit, and determining the status 
of these particles at the end of each run. There are either 
three or four possibilities, depending on the calculation: (a) 
the material has been involved in star formation; (b) the 


material has been ionised; (c) the material is neutral and 
denser than pcrit; (d) the material is neutral and less dense 
than Pcrit. In Figure [3] we plot bar charts for showing the 
results of this analysis, with the four categories of gas being 
shown as black, red, blue and green respectively. 

Runs I and J are globally bound clouds whose radii and 
therefore mean densities change little when undisturbed by 
feedback. In the control simulations of these clouds, the ma¬ 
jority of gas whose density at some point in the simulations 
comes to exceed pcrit is involved in star formation by the 
ends of the simulations. A small amount is still in gaseous 
form and denser than pcrit, and an even smaller amount 
has fallen to lower densities. In the feedback runs, it is 
evident that more dense gas is produced in total, but that 
less of it is eventually involved in star formation. However, 
the difference is not due to large quantities of dense gas 
being subsequently ionised, or to gas remaining at high 
densities but not forming stars. Instead, large quantities of 
gas are elevated over the density threshold, but later fall 
back below it, while remaining neutral. 

Runs UF, UP and UQ, are globally unbound and thus 
expand substantially whether feedback is acting or not. 
Consequently in the control simulations of these clouds, 
a substantial fraction (the majority in Runs UP and UQ 
in fact) of the gas whose density is raised above pcrit by 
turbulence (or gravity) subsequently re-expands and is 
never involved in star formation. We reiterate that this 
is nothing to do with feedback and is just a result of the 
initial supervirial states of these clouds. In the simulations 
in which feedback is active, the fraction of dense gas 
involved in star formation again drops and the total quan¬ 
tities of dense gas produced again increase, although only 
marginally in the case of Run UP. The fraction of dense gas 
which is ionised is once more very modest, although the 
ionised gas is the chief difference between the control and 
dual-feedback Runs UP. In Runs UF and UQ, feedback 
results in more dense gas being produced but much of 
it falls back below the density threshold without being 
ionised. 

These plots support several general conclusions: 

(i) The total quantity of sometime-dense gas in 
dual-feedback simulations is larger than in the control 
simulations. Feedback aids these clouds in generating dense 
gas. 

(ii) The fractions of dense gas which finishes simula¬ 
tions with densities lower than pcrit is always greater in 
the feedback calculations than in the corresponding control 
run. 

(iii) The quantities of dense gas which are ionised in the 
dual-feedback calculations are in all cases small fractions 
of the totals. Dense gas in these simulations is therefore 
mostly not prevented from forming stars by being ionised. 

(iv) The final star formation efficiencies are always 
lower in the dual-feedback simulations. Feedback restrains 
the clouds from forming stars. 

The obvious question arising from these observations 
is, if the feedback-affected clouds are better at forming 
dense gas, why are they worse at making stars? 
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(a) Control, tion (b) Control, tfin^l (c) Dual-feedback, tfl^ai 

Figure 4. The potential wells of the control simulation at tion (left) and tfinai (centre) and the dual-feedback simulation at tflnal (right) 
in Run I. Black dots are particles that are unbound in the cluster centre-of-mass frame and red dots are unbound particles. Bound 
particles with densities above Pcrit are shown in blue, and unbound particles exceeding this density are shown in green. Orange and 
teal triangles represent bound and unbound sink particles respectively. Unbound, in all cases, means having positive total energy in the 
system centre—of—mass frame. Note that only one tenth of all particles in each class are plotted. 



Figure 5. Gravitational potential against density at tjon (left) and tfjog^i (centre) in the control run, and in the dual-feedback simulation 
3-t tfinai ('"Ight) in Run I. Black dots are particles that are unbound in the cluster centre—of—mass frame and red dots are unbound 
particles. Bound particles with densities above pait are shown in blue, and unbound particles exceeding this density are shown in green. 
Unbound, in all cases, means having positive total energy in the system centre—of—mass frame. The vertical teal line is Pcrit- Note that 
only one tenth of all particles in each class are plotted. 


3.3 Conversion of dense gas to stars 

The answer to the question of why the feedback calculations 
are not efficient at converting their gas to stats lies in the 
fact that dense gas is a necessary but not sujficient condition 
for the production of stars. The gas must also be gravita¬ 
tionally bound and collapsing. In Figured we show the po¬ 
tential wells from the pre-feedback (left) and final (centre) 
states of the control Run I simulations, together with the fi¬ 
nal state of the dual-feedback simulation (right). Equivalent 
plots from runs J, UF, UP and UQ are very similar and are 
not shown to save space. Gas particles are shown as small 
dots. Only non-ionised gas particles are shown and, for rea¬ 
sons of clarity, only every tenth particle is plotted. Bound 
particles whose densities exceed pcrit are shown in blue, un¬ 
bound particles exceeding this density in green, and other 
unbound and bound particles as black or red respectively. 
Sink particles are shown as triangles, with orange indicat¬ 
ing bound sinks and teal unbound sinks. Several points are 
immediately apparent. 

(i) The potential well in the control simulation becomes 
deeper and broader with time as more gas falls into it and 
more stars are formed. The stars generally congregate in the 


deepest part of the well. In the final states of the control and 
dual-feedback simulations, there are a few unbound sink 
particles deep within the respective potential wells. These 
are low-mass sinks that have become unbound by dynami¬ 
cal interaction with more massive partners. 

(ii) Unbound gas particles in the control simulation are 
rather uniformly-distributed at values of the logarithm of 
the potential (expressed in kilometres per second) close to 
zero. In the ionised run, by contrast, unbound gas particles 
exist at substantially greater depths in the potential well, 
due to local acceleration of gas to velocities higher than the 
cloud escape velocity. 

(iii) Gas whose density is larger than pcrit resides almost 
exclusively in the deepest troughs in the potentials. This is 
particularly obvious in the control run. 

(iv) In the control simulation, the densest gas is coinci¬ 
dent with the densest groups of stars in the deepest potential 
troughs. However, in the dual feedback simulation, the deep¬ 
est potential well is devoid of gas and contains only stars. 
This well is substantially shallower than its counterpart in 
the final state of the control simulation. The dual-feedback 
simulation instead possesses several smaller shallower poten¬ 
tial wells containing mixtures of gas and stars, and not all 
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of the material in these wells is bound in the cloud centre- 
of-mass frame. 

Our interpretation of these observations is as follows. 
The cloud’s primary potential well in the control run is the 
main engine for the production of dense gas and for the 
conversion of gas to stars. The depth and dominance of this 
well ensures that gas continually flows into it and that this 
gas is bound and has little option but to form stars. In the 
feedback run, the reversal of the accretion flows feeding the 
central cluster and the expansion of the cluster itself has 
prevented the potential well acquiring the depth it achieves 
in the control simulation. Additionally, the deepest potential 
well in the feedback simulation has been completely emptied 
of gas, so it is not able to contribute to star formation. Com¬ 
pression of gas on the peripheries of the cleared-out bubbles 
has instead produced several less massive and shallower po¬ 
tential troughs. There is thus no large and efficient central 
star factory in the dual-feedback simulation and none of the 
gas has access to a potential as deep as that in the control 
simulations. 

The expansion of the HII regions/wind bubbles into 
the outer regions of the clouds in the dual-feedback sim¬ 
ulations does create additional gas whose densities exceed 
Pcrit for some time, but this gas never becomes gravitation¬ 
ally bound. We illustrate this in Figure [5] which depicts the 
potential of the gas as a function of density for the control 
run at the time of the initiation of feedback (left panel) and 
the end of the simulation (centre panel), compared to the 
end of the dual-feedback run (right panel). Black and red 
particles are unbonnd and bound particles respectively with 
densities below pcrit, while green and blue particles are un¬ 
bound and bound particles with densities above pcrit. 

In the control simulations at both epochs, all gas above 
Pcrit is bound in the centre-of-mass frame, and thus highly 
likely to form stars. In the dual-feedback run, by contrast, 
there are large quantities of unbound gas whose density ex¬ 
ceeds Pcrit. Although this gas is unbound in the frame of the 
cloud centre of mass, this does not necessarily imply that 
none of i t will become loc ally self-bound. However, as dis¬ 
cussed in lPale et al.l (l2013l l. because of the steep mean radial 
density profiles of the clouds, the surface density of 

swept-up material declines with time and radius and this 
gas becomes more gravitationally stable as the simulation 
progresses. We find that large quantities of gas which be¬ 
come denser than pcrit fail to become bound and eventually 
fall back below this density as they are driven to larger radii. 
This argues against a strict density threshold for star for¬ 
mation in these calculations. 


3.4 Stellar and star formation rate surface 
densities 

In iDale et ^ (120141 ') and in the preceding subsections, we 
discussed chiefly the global properties of the model clouds, 
such as their final star formation efficiencies and mean star 
formation rates. Here we examine the relations between 
resolved quantities, namely the stellar mass surface density 
or star formation rate surface density, and the gas mass 
surf ace density. _ 

iGutermuth et ahl (l201lf ) used resolved observations of 
eight nearby molecular clouds to examine the correlation 



Figure 6. Stellar mass surface density plotted against gas mass 
surface density for the final state of the control Run I simulation, 
at the location of each sink particle. Four points in a randomly- 
chosen colour are shown for each sink particle. In all cases, the 
stellar mass surface density is computed using the method of 
Gutermuth et al. 2011, with N=10. For the star symbols, the gas 
mass surface density has been computed by evaluating an SPH 
surface density sum at the locations of each sink particle. For the 
triangle, circle and square symbols, a uniform grid of pixels has 
been placed on a 15xl5pc box and the surface density computed 
by performing an SPH surface density sum at location of the pixel 
centres. Each sink particle has then been assigned the surface 
density of the pixel in which it s projected. Triangles represent 
results for a 64^ grid, circles a 128^ grid and squares a 256^ grid- 


between the stellar mass surface density and gas mass 
surface density. Gas surface densities are computed from 
near-IR dust extinction mapping, and the cloud total 
masses are computed from summing the mass at extinctions 
above detection thresholds peculiar to each cloud. The most 
massive system examined was Orion at 33 200 M©, and the 
least massive was Serpens at 2 590 Mq, with a mean of «14 
500 Mq, comparable to our low mass runs I, J, UP, UQ 
and UF. The Av thresholds used for these determinations 
varied from -1.0 (for Orion) to 3.7 (for Serpens), with a 
mean of « 0.7, or « 0.3 if Serpens is excluded. These figures 
illustrate the p roblem of comput i ng th e mass of a GMC 
observationally. IGutermuth et aP (1201 ill computed stellar 
surface densities by drawing a circle of radius i?N around 
each star cutting through the position of the A’th nearest 
neighbouring star, and computing the surface density as 
{N — 1)/'kR^, using N = 10. 

We replicate these techniques as closely as possible. We 
use the same method as iGutermuth et al.l ||201 ih to estimate 
the stellar density at the sink particle locations. We also 
note that typically 80-90% of the mass of our model clouds 
exists at Ay >0.5, so that their actual masses are close to 
what would be measured by extinction mapping with the 
typical thresholds used bv IGutermuth et al.l (1201 il l. 

However, it is not immediately obvious how to best 
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(c) 


(a) Run I 


i- - >P'^.' 

(b) Run UF 


Run E 


Figure 7. Stellar mass surface densities versus gas mass surface densities at the locations of all sink particles for Runs I (IO'^Mq bound 
cloud, left panel), UF (SxIO'^Mq unbound cloud, centre panel) and E (10® Mq bound cloud, right panel). Red triangles are from times 
when feedback is initiated, blue circles are from the end of the control runs, black squares are from the end of the dual-feedback runs. 
Hatched areas have A^f < 1, where observational results are likely to become unreliable. The polygons are from Gutermuth et al. 2009’s 
results; orange is Ophiuchus, purple is Mon R2, teal is Ceph OB2 and red is Orion. 



Figure 8. Star formation rate surface densities versus gas mass surface densities in Runs I (IO^Mq bound cloud, left panel), UF 
(SxIO'^Mq unbound cloud, centre panel) and E (10 ®Mq bound cloud, right panel). Pink triangles are from times when feedback is 
initiated, blue circles are from the end of the control runs, black squares are from the end of the dual-feedback runs. The black line 
shows the relation form Wu et al., 2005, the blue line is the Kennicutt et al., 1998 relation and the orange line is from Bigiel et al., 2008. 
The red region delineates approximately the region covered by the c2d and Gould Belt clouds as reported by Heiderman et al. (see their 
Figure 3). 


compute the gas mass surface densities, so we employ 
two techniques which are compared in Figure |6l using the 
control Run I simulation as a test. The results plotted in 
this figure are from the 2 -axis projection of the simulation, 
but we confirmed that observing along the x- or y-axes 
did not produce substantially different results. We first 
measured the gas surface density at the position of every 
sink particle by integrating through the smoothing kernels 
of every gas particle which overlaps in projection the 
position of the sink, to compute the projected mass as 
viewed along the 2 -axis. The same technique is used over a 
uniform grid to generate the column density images in our 
other papers. This technique is in some sense scale-free, 
since no grid is imposed on the density field. Alternatively, 
we placed uniform grids of 64 , 128 and 256 on the 


clouds and computed the column-density in each pixel by 
performing an SPH column-density sum at the location of 
each pixel centre. We then assign column-densities to the 
sink particles using the value of the pixel in which the sink 
is projected to lie. Figure [6] shows that these four different 
estimations of the gas surface density agree tolerably well 
over the range of grid resolutions tried. The 64^ grid is not 
able to capture the highest column densities, greater than 
~ 3 X IO^Mqpc”^, but otherwise the forms of the plots 
are rather similar. Since we are only going to use them 
for making qualitative judgements, we consider either of 
the techniques used adequate, and we adopt the first one, 
in which the column densities are computed at the sink 
particle positions. 

In Figure [T] we plot stellar surface densities against 
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gas surface densities for Runs I, UF and E, comparing the 
results at the epoch when feedback is enabled (red triangles) 
to those from the ends of the control runs (blue circles) and 
the ends of the dual-feedback runs (black squares). 

In all three pairs of simulations, the stellar surface 
densities grow in the control simulations relative to the 
epochs when feedback is enabled, by almost hve orders 
of magnitude in the case of Run UF. The maximum gas 
surface densities, by contrast, are generally somewhat 
lower at the ends of the control simulations, reflecting 
consumption of gas. The dual-feedback simulations all 
exhibit lower maximum stellar surface densities, due to the 
general result of feedback that the stellar systems have 
lower volume- and surface-densities, due to weakening of 
the local potential and decreases in the star formation effi¬ 
ciencies. The maximum gas surface densities are somewhat 
higher in the dual feedback simulations and tend not to 
be associated with high stellar surface densities, reflecting 
the relative inability of feedba ck to trigger star formation 
discussed in I Dale et al.l ll2013ah . 

W e compare our results to those of iGutermuth et aP 
ll201lli (their Figure 9). Included in Figure [7] are coloured 
polygons delineating approximately the areas of parameter 
space occupied by Ophiuchus (orange), Mon R2 (purple), 
Cepheus OB2 (teal) and Orion (red). The control simula¬ 
tions, particularly of Runs I and UF, resemble most in form 
Ophiuchus and Mon R2, in that generally higher stellar 
surface densities correspond to higher gas densities. This 
correlation is quite tight in the case of the observations, 
as illustrated by the orange and purple polygons, and the 
relation between stellar mass and gas mass surface density 
follows approximately E* oc and E* oc Eg^®^. 

The control Run I simulation initially exhibits a tight 
correlation between E* and Egas but the relations at the 
ends of the control runs all show more substructure and 
considerably weaker correlations. The formal fitted slopes 
for the final states of the control runs I, UF and E are 
0.99, 2.58 and 1.70 respectively, so that the slopes of the 
E,/Egas relations for runs UF and E are similar to those 
of Ophiuchus and Mon R2. However, we also observe sub¬ 
stantially higher maximum stellar mass s urface densities in 
the co ntrol simulations than are shown bv iGutermuth et al.l 
(l201lli . This likely has several causes. Resolution limitations 
in the observations underestimate the surface densities of 
the densest regions in the real clouds, but Ophiuchus, being 
a low-mass star-forming region is not expected to exhibit 
very high stellar surface densities in any case. However, 
probably more importantly, our control simulations where 
star formation is entirely unrestrained produce too many 
stars in configurations which are very compact, leading to 
very high stellar surface densities 

The feedback simulations, in contrast, produce sub¬ 
stantially lower s tellar surface den s ities more similar to 
those observed by IGutermuth et al.l (l201ll l. particularly in 
the case of Runs I and UF, where the maximum stellar 
densities are reduced by approximately two orders of 
magnitude relative to the control runs at the same epoch. 
The f eedback calculations resemble more IGutermuth et aH 
(l201lli ’s Cepheus OB3 and Orion data, illustrated by the 
teal and red polygons respectively. The simulations and 
observations exhibit large horizontal spreads, correspond¬ 
ing to regions with high stellar densities that have been 


partially or largely cleared of gas. Termination of accretion 
flows, stifling of cluster growth, and spatial separation of 
stars and gas in the feedback simulations produce stellar 
surface densities and correlations between stellar and gas 
mass surface densities closer to those from Gutermuth et 
al. i n systems where feedba ck active. 

iHeiderman et all (I2OIC1II instead plotted the star forma¬ 
tion rate surface density for twenty nearby c2d and Gould’s 
Belt clouds with masses ranging from 189 Mq (Lupus IV) 
to 24 400 Mq (Serpens-Aquila). They also use extinction 
mapping, complemented with ^^CO and observations. 

They divide the clouds into regions enclosed between 
extinction contours and compute gas and star formation 
rate surface densities for each contour for each cloud. They 
compute the star formation rate surface densities from YSO 
surface densities, assuming a mean YSO mass and age. 

We also partition our clouds into contours based on 
extinction. We choose a minimum Av of 4 magnitudes 
and a contour spacing of 4 magnitudes, with a maximum 
number of contours of 8. Within each region defined by two 
contours, we compute the mean gas surface density. Since 
our sink particles are not all of the same mass and we have 
time information available, we obtain star formation rate 
surface densities in each region by locating all sink particles 
within that region and referring back to previous dumps 
0.5Myr in the past to compute the total gain in stellar mass 
in each region. 

We plot the results in Figure [S] again for Runs I, 
UF and E, with pink triangles from times when feedback 
is initiated, blue circles from the end of the control runs 
and bla ck squares from t h e end of the dual-feedback runs. 
As did IHeiderman et af] (I 2 OIC 1 IL we also include the star 
formation rate surface de n sity v er sus gas su rf ace d ensity 
relations from Wu et ahl (l2005h : iKennicut^ lll998h and 
iBigiel et al.l ll20 08h to mide the eye. No t e, how ever, that 
the iKennicuttl 1 199811 and iBigiel et al.l ll2008h relations 
are derived at much larger size scales and thus effectively 
smear out star formation activity over large areas. They 
are thus not expected to give an accurate representation of 
star formation at the level of indivi d ual m olecular clouds. 
The results from IHeiderman et all ll2010ll . however, are 
derived at mass and scales similar to our simulated clouds 
or to subregions of them, so we plot red polygons which 
approximately delineate the area occupied by their c2d and 

Gould’s Belt results (see their Figure 3). _ 

_ Co mparing our res ults t o the mode l s of IWu et al.l 

1I2OO5II ; iKennicu^ lll 998 ll and IBigiel et al.l 1I2OO8II , we Hnd 
that our star formation rate densities are almost always 


substantia ll y high er a t a given den si ty tha n predicted by the 
iKennicut^ Il99^ or IBigiel et al.1 (I 2 OO 8 II relations. Again, 
this is likely due in part to the fact that these relations are 
derived from observations at scales of ~ lOOpc and above, 
so that star formation is not well resolved. At g as surface 
densit ies close to IO^Mqpc”^, the agreement with lWu et al.l 
ll2005h is rather better, but at higher gas surface densities, 
we again recover much larger star formation rate surface 
densities in Runs I and E, although the agreement in Run 
UF is rather closer. 

Our results are in fac t much closer in form to those 
of IHeiderman et al.l ll2010ll ’s Figure 3, exhibiting a much 
steeper relationship between Esfr and Egas. Runs I and 
E have substantially higher star formation rate surface 
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densities at a Riven ga s surface density than reported by 
iHeiderman et all ll2010l '). but Run UF is rather closer. This 
is likely due to Run UF being an unbound cloud with, on 
average, lower star formation rates than the two bound 
clouds. 

In all three simulations (which are representative of 
the others), we see that the differences between the control 
simulation at either tion or tfinai and the dual-feedback 
simulations at tfinai are modest. All three sets of points 
occupy broadly the same region in the plots and there are no 
clear features to distinguish the control and dual-feedback 
calculations. It is thus difficult to discern the influence of 
feedback simply by analysing correlations between star 
formation rate surface density and gas surface density in 
this manner. This is in contrast to the results presented 
earlier where there is a clear effect of feedback on the 
correlation between stellar surface density and gas surface 
density. Star formation occurs in dense gas regardless of 
whether feedback is active or not but, as we saw in earlier 
sections, not all the dense gas in the feedback simulations 
is active in forming stars. Dense gas which is not forming 
stars is obviously going to be poorly represented in plots 
constructed in the manner of Figure[8l based on YSO counts. 


3.5 Unbinding of stars and clusters 

In Figure [9] we show the gas mass fraction and mass 
and number fractions of stars unbound at the ends of 
our feedback calculations as a function of cloud escape 
velocities, expressed in units of the sound speed inside the 
HII regions (11 km s“^). The relevant unbound fractions 
from the control simulations have been subtracted to isolate 
the effects of feedback, so the plots effectively show the 
quantities of material unbound by the action of feedback. In 
the case of clouds where the mass or number fractions of 
stars unbound is very small, the mass or number fractions 
unbound in the control simulations may be slightly greater 
than in the feedback simulation. Therefore, not all simula¬ 
tions appear in the centre and right panels of the hgure. 

There is a tight correlation between the unbound 
gas mass a nd the cloud escape velocity, as discussed in 
IPale et ^ (l2013bh . The correlation between the unbound 
stellar numbers and the cloud escape velocity is weaker 
with more scatter, and several clouds absent, from having 
greater numbers of stars unbound in the control simulations 
that in the dual-feedback run. There appears to be no 
correlation at all between the unbound stellar mass and the 
cloud escape velocity. We check this in Figure [TOl where we 
show that the unbound stellar mass fraction is apparently 
uncorrelated with but generally lower than, the unbound 
gas mass fraction, so that combined photoionisation and 
wind feedback is generally much less effective at unbinding 
stars than it is in unbinding gas. We explore the reasons for 
this lack of correlation in the Discussion section. 

Figure [9] suggests that neither the disruption of the 
clouds nor that of the clusters is mass-independent. In 
the case of the clouds , one im portant cause of this, as 
discussed in IPale et ^ ll2013bl l. is likely to be the cloud 
escape velocity, or binding energy. In order to unbind the 
clouds, one must either supply sufficient momentum pcrit 
to accelerate the whole cloud beyond its escape velocity, or 



10“^ lo-* 

Unbound gas mass fraction 


Figure 10. Unbound stellar mass fractions plotted against un¬ 
bound gas mass fractions. Colours again denote cloud masses, 
circles are for initially unbound clouds and squares are initially 
bound clouds. In all cases, the relevant fractions from the control 
simulations have been subtracted, so not all simulations appear 
in this plot. 


energy £?crit in excess of its binding energy. Since the clouds 
all have roughly the same surface densities, M so 

VESC - Pcrit ~ and Scrit ~ However, 

as also mentioned in our earlier paper, the ability of the 
feedback-driven bubbles to explore the cloud volume also 
strongly affects their ability to destroy the clouds. At the 
ionised sound speed, an HII region can expand at most 
« 30pc in the 3 Myr time window we have allowed here 
before the detonation of supernovae, which is greater than 
or comparable to the radii of the lower mass clouds, but 
much smaller than the radii of the more massive clouds. 

The vulnerability of clouds to feedback s hould also 
depen d on their stellar content. As discussed in iDale et al.l 
ll2014ll and shown in the leftmost panel of Figure 1111 with 
the exception of the nearly gas-exhausted Run F, the final 
star formation efficiencies of the simulations span about one 
decade from 1.3 to 11%. Computing for each simulation the 
fraction of the cloud mass in stars above 20Mq (by simply 
counting in the 10"^ and 3xlO‘^M0 clouds and by assuming 
each subcluster in the more massive clouds has a Salpeter 
mass function between 0.5 and 100 Mq and integrating for 
the more massive clouds) reveals a similar spread across 
the parameter space (centre-left panel of Figure [TTJ. 

From the point of view of feedback, a more important 
quantity is the ionising luminosity per unit mass of the 
clouds (since ionisation dominates over winds), shown in the 
centre-right panel of Figure 1111 which hints that the more 
massive clouds have smaller ionising luminosities per unit 
mass, although again with large dispersions, which is partly 
accounted for by the unbound clouds having generally lower 
star formation efficiencies. However, it may also be due to 
the simplistic method used to assign feedback luminosities 
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Figure 9. Unbound gas mass fractions (left panel), stellar number fractions (centre panel) and stellar mass fractions (right panel) plotted 
against cloud escape velocities normalised to the ionised sound speed (11 km s“^) for all simulations. Colours denote cloud masses, circles 
are for initially unbound clouds and squares are initially bound clouds. In all cases, the relevant fractions from the control simulations 
have been subtracted, so not all simulations appear on all plots. 



Simulation mass (My) Simulation mass (My) Simulation mass (M „) Simulation mass (M.„) 


(a) (b) (c) (d) 

Figure 11. Final star formation efficiencies (leftmost panel), fraction of simulation mass in O-stars (centre-left panel), numbers of ionising 
photons per IO^Mq (centre—right panel) and final global ionisation fraction (rightmost panel) for the final state of all simulations, with 
initially bound clouds shown as teal squares and initially unbound clouds as purple circles 


to the clouds of masses IO^Mq and above necessitated 
by our inability to resolve individual stars in these runs. 
This underestimates the luminosities of the most massive 
subclusters because it treats 30 Mq stars as the basic unit 
of feedback, ignoring possible contributions from rare but 

much brighter more massive o bjects. _ 

However, we showed in iDale et ^ ||2012| ) that the 
results of simulations are not sensitive to changes of factors 
of several in the ionising luminosity. As we show in the 
rightmost panel of Figure ITTl the dispersion in the fraction 
of clouds which are ionised is quite small with all calcu¬ 
lations (save Run F) having between 1 and 12 % of their 
mass io nised. Inspired by this, we presented in iDale et ^ 
ll2013bl ) a simple model treating the HII regions as pistons 
and considering what fraction of the cloud mass could be 
raised by their expansion to the appropriate escape velocity, 
finding a strong dependence of this fraction on the escape 
velocity itself. 


3.6 Cluster emergence 

The observed emergence of stellar clusters from their em¬ 
bedded phase occurs as a result of gas being cleared from 
the clusters, regardless of whether the loss of gas dynami¬ 
cally impacts the clusters. A detailed numerical study of this 
phenomenon strictly requires advanced radiation transport 
calculations which are beyond the scope of this work. How¬ 
ever, some insight can be gained by a simple analysis. 

Open clusters and stellar associations are generally eas¬ 
ily visible in the optical (e.g. the V-band), whereas embed¬ 
ded clusters typically lie behind 5-100 ma gnitudes of extinc¬ 
tion in the V-band llLada fc Ladal [20031) and are rendered 
nearly invisible at optical wavelengths. It is the transition 
between these states which is generally referred to as ‘emer¬ 
gence’. Infrared (e.g. K-band) extinctions are around eight 
times smaller than optical extinctions for a given column 
density, so embedded clusters obscured in the optical may 
be easily visible in the infrared. 

Protostars and pre-main-sequence objects are very 
bright in the infrared thanks to their accretion luminosi¬ 
ties and the reprocessing of radiation by their dust-laden 
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gaseous envelopes. Modelling the K-band luminosities of 
such objects is very difficult and we will not attempt it here. 
The visible and ultraviolet light, however, should be domi¬ 
nated by the bright massive stars, which become well settled 
on the main sequence while still embedded. We therefore 
choose to treat the massive stars in our clouds as having 
main sequence colours and luminosities and examine the ex¬ 
tinction of their UV/visible light by the remaining gas in the 
simulations. We then sum the contribution of the massive 
stars to estimate the total visible magnitude of the whole 
stellar population, as it would appear if it were an unre¬ 
solved cluster. 

We utilise tabulated colour s computed for stars i n Hub - 
ble Space Telescope filters from iLeieune fc Schaereil ll200lli . 
We choose solar metallicity and IMyr isochrones as a reason¬ 
able average age for our stars. We concentrate on the F336W 
(approximating the U-band), F439W (the B-band) and 
F555W (the V-band) Hlters. We include contributions from 
all stars more massive than lOM©. We compute the intrinsic 
magnitude of each star from its mass only and obtain the 
extinction in front of each star from the column density be¬ 
tween the star and an observer at infinity in the x-, y- or z- 
directions. The column densities are converted to extinctions 
by assuming that U-, B- and V-band extinction laws apply, 
so that Av=VH/10^^cm“®, Ab = 1.322Av (this correspond s 
to the canonical mean extinction law fe.g. lFitzpatricklll999ll . 
but we note that values as low as AB~l-322Ay may be ap- 
*ropriate for dense reg ions of the ISM) and ^u=1.581j4v 
S avage fc Mathi^ligTS l. 


Figure [T2] shows, for the final timesteps of the control 
(blue) and dual-feedback (grey) Run I simulations, the val¬ 
ues of Ay for all stars when viewing the system along the 
a-axis, plotted as a function of stellar mass. The extinctions 
in the control run are sharply peaked at Ay ~ lOmag and 
extend up to «50 mag. Those in the dual-feedback run fall 
into two clear groups. One, comprising roughly half the stars 
and corresponding to the population of stars embedded in or 
near the bubble walls, has a mean extinction also close to 10 
magnitudes. The other half have extinctions less then 1 mag 
and often negligible, corresponding to the exposed cluster 
inside the cleared out feedback-driven bubble. Apart from 
the few very most massive stars in the feedback simulation, 
whose extinctions are very low, there is little dependence of 
extinction on stellar mass in either simulation. In particu¬ 
lar, in the feedback simulation, a star of any mass can have 
its environs completely evacuated of gas by the few massive 
stars present in the system. 

To compute the magnitudes of the clusters, we first 
convert the magnitudes of each star to fluxes in arbitrary 
units using 


F'c1us,u,b.v = XI E , v , B,y = X ^0 ^ io-(^u,b,v-Au,b,v)/2.5^ 

stars stars 

(5) 

then sum over all sources and convert back to magnitudes 
using 

M<,iu,,u.b,v = -2.51og ^ (g) 

In Figure we plot the total magnitude in the F555W 
Hlter as a function of time for the control (blue lines) and 
dual-feedback (black lines) Runs I (left), J (centre) and UQ 
(right) simulations. In each case, we plot results looking 
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Figure 12. Foreground extinctions in V-magnitudes for all stars 
in the control (blue) and dual-feedback (grey) as a function of 
stellar mass in Run I at the ends of the simulations, at times of 
7.56Myr when viewing the system along the 5:—axis. 


down the x- (dash-dot lines), y- (dashed lines) and a-axis 
(solid lines) to gain an idea of how projection effects may 
influence the results. 

The iuitial magnitudes for the three clusters are in the 
range 5-25. There is substantial variation in the results de¬ 
pending on the viewing angle. In Run J in particular, chang¬ 
ing the viewing angle results in variations of almost 20 mag¬ 
nitudes in the cluster brightness. This is purely due to the 
highly anisotropic nature of the dense gas in the which the 
majority of the stars are initially embedded. However, the 
same qualitative conclusions may be drawn from all three 
viewing angles. 

In both the control and dual-feedback simulations, the 
clusters generally become brighter with time, often quite 
rapidly. In the control simulations, this is due largely to un¬ 
abated accretion producing brighter and brighter stars. The 
increase in brightness is not smooth, however, with consid¬ 
erable jumps visible. These features are a result of the non¬ 
steady delivery of gas to the core regions of the cloud by ac¬ 
cretion flows along the filaments, leading to large excursions 
in extinction. In the dual-feedback simulations, by contrast, 
the increase in cluster brightness is generally much more 
abrupt, with all three of these simulations producing ex¬ 
posed clusters with brightnesses of -6- -7 mag in the F555W 
filter, with differences in viewing angle resulting in spreads 
of about 1 mag. However, the timescales on which the clus¬ 
ters are revealed also depend somewhat on viewing angle. 
In Run J, the cluster as viewed along the z-axis appears to 
achieve a magnitude brighter than zero only «10®yr after 
ionisation is enabled, whereas this timescale is R;4xl0®yr 
for an observer on the y-axis. The greater variation of this 
timescale in Run J is due to the less effective clearing out 
of the cluster volume by feedback in this calculation relative 
to the other two shown. 
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Figure 13. Absolute V-band magnitudes of clusters as a function of time in Runs I (left panel), J (centre panel) and UQ (right panel). 
Control simulations are shown as blue lines and dual-feedback runs as black lines. Each cluster is observed along the x— (dash-dot lines), 
y- (dashed lines) and 2 :-axis (solid lines). 


Run 

tSF 

tpB 

t-5 

^end 
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ffb 
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I 
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0 39+0 07 
u.oy_Q ig 

0.59l°i2 

UF 

2.08 

3.28 

4 43+0.13 

6.23 

4.15 

2.95 
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n 43+^-^*^ 

Q3 

U.D1_0.04 

UP 

0.92 

1.83 

z.du_o.72 

3.71 

2.69 

1.88 

l-lllale 

0 41+01’^ 

u.oy_o.24 

UQ 

1.94 

3.13 

Q >70+0.99 

O- '0-0.64 

5.86 

3.92 

2.73 


n 53+016 
u.oo_Q 26 

n 7fl+6-24 

U. ' 0 _Q 37 


Table 2. Comparison of the times (all in Myr) at which star formation begins (tsF)i the time at which feedback is enabled (tFB)i the 
time at which the clusters reach V=-5 (t— 5 ), the end time of each simulation (tend)i the total duration of star formation tsfi the time 
since feedback was enabled tfbi the time for which the cluster in the feedback run has been brighter than a V—magnitude of -5 t_ 5 , the 
fraction of time since star formation began for which the clusters have been brighter than V=-5 (fsF) and the fraction of the time since 
feedback was enabled for which the clusters have been brighter than V=-5 (fps)) for Runs I, J, UF, UP, UQ. 


Figure [TJ] shows the evolution of the cluster colour- 
magnitude diagram in the F555W versus F439W-F555W 
space of the control (blue-green colour map and circles) and 
dual-feedback (yellow-red colour map and triangles) Run I, 
as viewed along the jz-axis, with the colour bar at bottom 
left giving the elapsed time since feedback was enabled in the 
dual-feedback calculation. The control simulation becomes 
brighter and bluer but, at later times, oscillates about a 
colour of « 1.5 and a magnitude of ~ 0. By contrast, the 
dual-feedback simulation moves very rapidly from the bot¬ 
tom right to the top left of the CMD, quickly acquiring a 
stable colour of « —0.3 and stable magnitude of « —7. 

Figure [T^ shows a CMD for the individual stars (> 
IOM 0 ) in the control and dual-feedback Run I simulations 
at the ends of the runs. All of the stars in the control simu¬ 
lation are strongly reddened and control run stars of a given 
colour are several magnitudes brighter than stars of the 
same colour from the dual-feedback simulation. The stars 
from the control simulation are intrinsically brighter, since 
they are able to acquire larger masses, but are buried in 
much larger extinctions, so appear redder. The population 
embedded in the bubble walls in the dual-feedback simula¬ 
tion is visible at a colour of 2-3 and a magnitude of « 5. 
The exposed population is visible at the top left as a ver¬ 
tical grouping characteristic of massive stars on the main 


sequence, whose colours vary little but whose magnitudes 
vary substantially. 

The timescale on which clusters become observable is of 
crucial observational importance and can be estimated using 
the above analysis. In Table O we give, for Runs I, J, UF, 
UP and UQ the time at which star formation began (Isf), 
the time at which feedback was enabled (tfb), the time at 
which the clusters became brighter than V=-5 (t-s), chosen 
in part because no cluster in any control simulation achieves 
this brightness) and the time of the end of the simulations 
(tend). We use these figures to compute, at the endpoints 
of the simulations, the total duration of star formation tsf, 
the time since feedback was enabled tfb , the time for which 
the cluster in the feedback run has been brighter than a V- 
magnitude of -5 t _5 , the fraction of the star-forming life¬ 
time for which the cluster is brighter than V=-5, fsp, and 
the fraction of the time for which the O-stars have been 
active for which the cluster is brighter than V=-5, fpB. It is 
clear from Figure[T4]that t _5 varies substantially with view¬ 
ing angle. We show this crudely by giving for each cloud the 
median value of t-s from the x-, y- and 2 -projections, with 
the other two as rough upper and lower limits. 

As alluded to previously, the variation with viewing an¬ 
gle in the times at which the clusters become brighter then 
V=-5 can be substantial, with variations up to almost IMyr 
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Figure 14. Evolution of the colour—magnitude diagram of the 
clusters from the control (blue—green colour map and circles) 
and dual-feedback (yellow-red colour map and triangles) Run 
I calculations in the F555W versus F439W-F555W space. The 
elapsed time since the point when feedback is enabled in the dual¬ 
feedback run is given in the colour bars at bottom left. 


in the case of Run UQ. However, the time delay between 
the onset of star formation and the formation of the first 
few O-stars, defining when feedback is enabled, is also ap¬ 
proximately IMyr in all simulations, and the star formation 
timescales are in the range 2-4 Myr in total (and would be 
longer in the case of Runs J and UP if it had been practi¬ 
cal to continue the simulations for the full 3Myr after the 
initiation of feedback). The fractions of the star formation 
timescales for which the clusters are brighter than V=-5 are 
in the range 39-54% with typical variations of « 15% and 
absolute maxima and minima 69% and 23% respectively. 
Clusters should thus be bright for about half the time in¬ 
terval between the onset of star formation and the detona¬ 
tion of the first supernovae, but the scatter in this quantity 
is significant. Naturally, the fraction of time since the O- 
stars formed for which the clusters are this bright is higher, 
between 59 and 84% with absolute maxima and minima of 
100% and 25%. The variations in this quantity are thus even 
larger. These variations imply that substantial fractions of 
cluster in magnitude-limited surveys may be missed purely 
through inhomogeneities in the intra-cluster gas. 


4 DISCUSSION 

4.1 Disruption of clouds and clusters 

We found that the disruption of the clouds by feedback 
is not mass-independent and that the disruption of the 
stellar clusters is also not mass-independent, and is sub¬ 
stantially less severe. T his result diffe r s from previous 
work on this issue (e.g iTutukovI Il978l : iBoilv fc Kroupal 


Figure 15. Colour-magnitude diagram of individual stars in the 
control (blue circles) and dual-feedback (black squares) Run I 
calculations in the F555W versus F439W-F555W space at the 
endpoints of the two calculations. 


2003bl: Goodwin & Bastian 

2006; 

20071: Ipfalzner & KaczmarekI 

2Q1(; 


iBaumgardt fc Krounal 
20131) . These authors model 


the effects of gas expulsion by implicitly or explicitly 
allowing an N-body system to come into equilibrium 
with an artificial potential, then removing that potential 
instantaneously or on a prescribed timescale. The stellar 
system modelled is initially smooth, often in the form of 
a Plummer sphere, and the background potential is also 
smooth. The contribution of the stellar mass to the total 
mass density (the ‘star formation efficiency’) is taken to be 
the same at all radii. 

These are tempting simplifications, but hydrodynamic 
simulations in which star formation is self -consistently 
model le d have shown th a t they are not correct. Offneret al.l 
(l2009l) : iKruiissen et al.l (12012) : iGirichidis et al.l ll2012l) all 
find that the stellar and gaseous components of embedded 
clusters are never in virial equilibrium with each other, 
and that the stars tend to be subvirial in both the absolute 
sense of having a virial ratio less than 0.5, and a relative 
sense in having a smaller virial ratio than that of the gas. 
The stars and gas thus form two largely decoupled systems. 

We find that the same is true in our simulations. Figure 
1161 depicts, for Run I, the virial ratio of all cold gas (solid 
lines), the cold gas occupying the same volume as the stars 
(defined by the smallest sphere centred on the most massive 
which contains all the stars, dotted lines), and of all stars 
(dashed lines) in the control (blue lines) and dual-feedback 
(black lines) calculations. 

We here define the virial ratio as a simple ratio of ki¬ 
netic energy to gravitational potential energy i5kin/|F^grav|. 
The kinetic energy for each component (e.g. cold gas, 
stars) is computed by summing the kinetic energies of each 
individual star or gas particle in the centre of mass frame 
of the whole simulation. The potential energy is computed 
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Figure 16. Virial ratios of all cold gas (solid lines), the gas oc¬ 
cupying the volume defined by the location of the most massive 
star and the furthest star from this point (doted lines), and of all 
stars (dashed lines) in the control (blue lines) and dual-feedback 
(black lines) Run I calculations. 



Distance from most massive star (pc) 


Figure 17. Fractional contribution of stellar mass to the total 
mass in radial bins as a function of radius from the most massive 
star at the onset of ionisation (red line and triangles), the end 
of the control run (blue line and circles) and at the end of the 
dual-feedback (black line and squares) in Run I. 


by summing over all star or gas particles the total potential 
energy generated by all mass components of the simulation 
on each particle. 

The stellar virial ratio is always less than that of the 
gas in both calculations. In the control simulation, the 
virial ratios of both components change only modestly 
over the course of the calculation, the stars becoming more 
subvirial, while the gas slowly approaches a virial ratio of 
unity, likely due to the slow expansion of the outer regions 
of the cloud. In the feedback simulation, the virial ratio of 
the cold gas increases strongly and rapidly, but that of the 
stars increases much more slowly. 

If only gas inside the volume occupied by the stars is 
considered, the picture is slightly more complex. Initially, 
the virial ratios of the stars and cospatial gas are very 
similar. When feedback is enabled, the virial ratio of this 
gas in the dual-feedback simulation rises very sharply, since 
it is quickly and directly affected by feedback. As time 
progresses, the virial ratio of this material initially drops, 
before rising again at late times as the second-generation 
feedback sources become active. However, at all times after 
feedback is enabled, this gas has a substantially higher 
virial ratio than the stars. In the control simulation, the 
virial ratio of the gas cospatial with the stars simply rises 
gradually with time. Since the stars become more subvirial, 
the intracluster cluster becomes somewhat supervirial with 
respect to the stars. 

The virial states of the gas and stars in these simula¬ 
tions are never strongly coupled, even when only considering 
the gas cospatial with the stars. This is especially true of 
the dual-feedback simulation, where gas inside the cluster 
is always supervirial with respect to the stars by a factor of 
at least four. It is not surprising that the expulsion of even 
a large fraction of the surviving gas in some of the feedback 
simulations has only a modest effect on the dynamics of the 
stars. 

We reinforce this point in Figure 1171 which shows the 
fractional contribution of stellar mass to the total mass in 
radial bins as a function of radius from the most massive 
star in the control run at the onset of ionisation (red), 
at the end of the control run (blue) and at the end of 
the dual-feedback Run I (black). The total star formation 
efficiencies averaged over the entire clouds are, respectively, 
4, 13 and 8 percent. However, the local star formation 
efficiencies clearly depart strongly from the average. This is 
particularly true in the dual-feedback simulation, where the 
clearing of gas from the central cluster and its expansion 
result in a region «6pc in radius where all the mass is 
stellar. However, even in the absence of feedback, the 
local star formation efficiencies in the much more compact 
central cluster both at the onset of ionisation and at the end 
of the control simulation are much higher than the global 
average. The assumption that the gas and stars are evenly 
mixed is therefore a poor one. The local star formation 
efficiency in the central cluster at the onset of ionisation is 
wel l in excess of the critical valu e of 33 percent identified 
by iBaumgardt fc Krounal ll2007l ) as being the minimum 
required to survive adiabatic gas expulsion, so the survival 
of this stellar system in the dual-feedback simulation is 
not, in fact, surprising. Although more than half of the 
cloud mass is unbound by feedback, the gas is mostly lost 
from regions where there are no stars. 
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Recently, more advanced N-body simulations have 
improved on the methodology of early cal culations. These 
enhan c ements have taken two main forms. iMoeckel fc Bat j 
ll2010ll : iMoeckel et al.l ll2012h dropped the assumption of 
smooth stellar distributions by extracting the sink p articles 
for med in the simulatio ns of (respectively) iBa 3 (|2009all 
and iBonnell et ahl (l201l[l . By dispensing with the gas, they 
effectively assumed instantaneous gas expulsion, which 
is the most des t ructiv e to the remaining stellar system. 
IMoeckel fc Bate _|20K ) foun d that 30-40 percent of the 
stars from the Bate! ( 2009al l calculation remained bound 
after 10 Myr of post-gas expulsion evolution, despite 
the hydrodynamical calculation having a star formation 
efficiency of 38 percent, well below the can onical 50 p ercent 
required t o survive instantaneou s gas loss (|Hillslll98nl'l . The 
results of IMoeckel et al.1 (|2012| '1 were qualitatively similar 
and fo r similar reasons - th e clusters formed in the calcula¬ 
tion of iBonnell et al.1 (1201 ill were already gas-d epleted and 
viri alised, as report e d by Kruiissen et al. J 2012ll. 

IParker fc Paid (12013 1 and Parker et al. ( 2015l l used a 
si milar appro a ch to c ontinue a s u bset o f the simulations 
of Dale et all (I 2 OI 2 I I. iDale et ^ (l2013bl l and iDale et ^ 
(I 2 OI 4 I 1 . They assumed that, after the «3Myr action of 
ionisation feedback from the clusters’ 0-stars, the first 
supernova in each model cloud instantaneously expelled the 
remaining gas. By comparison with control simulations in 
which no pre-supernova feedback operated, they examined 
what effect the slow gas expulsion caused by ionising 
feedback might have on the evolution of the clusters after 
sudden but delayed gas expulsion. They found that it 
was difficult to identify clear trends in the effects of pre¬ 
supernova feedback, but that in all cases, the instantaneous 
removal of the remaining (usually large) fraction of gas left 
substantial fractions of the stars bound (more than half in 
all but one of the ten simulations analysed). The fractions 
of unbound stars in all calculations increased with time 
over the 10 Myr of N-body evolution, but this was a result 
of stellar-dynamical and stellar-evolutionary effects. 

Other authors have sought to improve the N-body 
modelling of young clusters by using more r e alistic , but 
still artificia l , init ial conditions. ISmith et al.l (1201 il l and 
ISmith et al.l (I 2 OI 3 I I construct artificial substructured clus¬ 
ters using either a fractal space-filling method, or by 
placing a number of small Plummer spheres within a larger 
Plummer potential. This has the effect of ensuring that 
the local star formation efficiency can be much higher than 
the mean efficiency, as we observed in Figure [T71 and these 
authors find that the local SFE is the main determinant 
of cluster survival, provided that gas expulsion occurs 
late enough that the stars have relaxed. However, the 
background potential representing the gas is still smooth 
and can still only be homologously removed. 


4.2 Emergence of embedded clusters 

For the low-mass clouds, the stellar populations in the 
control simulations never become fully optically revealed, 
although they do become substantially brighter over the 
course of the simulations, due to unrestrained accretion. 
Their apparent brightness can vary by several magnitudes 
on timescales of order 10®yr, due to non-steady flows of gas 


into the central regions of the clouds. In the dual-feedback 
simulations, we find that the timescales on which the 
clusters emerge from their veiling gas and become optically 
visible depends on the structure of the gas in the clouds 
and locate of the stars within it, but is short, generally 
0.1-0.5Myr and never more than IMyr. This timescale is 
substantially shorter than that on which the clouds are 
being disrupted. The times in which clouds I, J and UQ 
have half of their mass unbound are >2Myr. Although 
many of the stars in these simulations remain in regions 
of high extinction, the most massive and brightest objects 
which dominate the visible luminosity are very quick to 
disperse their surrounding material and become visible. 
This implies that the timescale for cloud dispersal and 
for the optical emergence of embedded clusters are not 
necessarily strongly connected. 


5 CONCLUSIONS 

We have evaluated the combined effects of photoionisation 
and wind feedback from O-type stars on the young embed¬ 
ded clusters formed in model turbulent GMCs o f a range of 
mass es and radii. We showed in previous work JPale et al.l 
l2014h that the fraction of gas expelled from the clouds was 
strongly dependent on the cloud escape velocities, but here 
we concentrated on the effects of feedback on the clusters, 
the ongoing star formation process and the interplay 
between gas and stars. Our conclusions are outlined below: 

(i) In the control simulations without feedback, star 
formation proceeds at the freefall rate in the dense gas, so 
that the depletion time of this material corresponds to its 
freefall time. In the feedback simulations the star formation 
rate becomes decoupled from the freefall time in the dense 
gas, sometimes being over an order of magnitude slower. 

(ii) The manner in which the star-formation- 
rate/dense-gas-mass relation is broken does not involve 
very large changes in the star formation rate, but instead 
is due to the production of dense but non-star-forming 
gas by feedback. Feedback expels dense gas from the 
clouds’ potential wells, and generates dense gas outside the 
potential wells which is less able to form stars, but only 
affects the star formation rates by factors of at most two. 
If regions where star formation is occurring within clouds 
are spatially correlated with dense gas, the relationship 
still appears to be strong, as shown in Figure [H since star 
formation always takes place in dense gas. Instead, feedback 
creates large quantities of dense gas which fail to form stars. 

(hi) The increase in dense gas but the similar or lower 
star formation rates in dual-feedback simulations relative 
to control simulations is due to the fact that most dense gas 
in the feedback calculations resides in the outskirts of the 
clouds where expanding bnbbles have swept up low-density 
quiescent gas. The dense gas in these rnns is not to be found 
in the depths of the clouds’ potential wells which are, in the 
control simulations, the principal star-formation engines, 
serving to concentrate the dense gas further and preventing 
it from dispersing. Secondarily, the expulsion of gas from 
the potential wells and suppression of accretion flows 
feeding the clusters in the feedback simulations reduces the 
depths of the potential wells relative to those in the control 
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calculations. 

(iv) Correlating star formation rates with gas surface 
density does not reveal strong differences between the 
feedback and control simulations. These quantities show a 
strong relationship in both cases. Conversely, correlating 
stellar surface density and gas surface density does reveal 
such differences. The separation of stars and dense gas 
caused by the clearing out of the central clusters and the 
generation of non-star-forming dense gas in the feedback 
simulations produces markedly different behaviour from 
the control simulations, where stars and dense gas remain 
correlated. 

(v) Simple modelling of the observed emergence of the 
clusters in the near-UV and optical bands implies that the 
apparent magnitudes of the clusters depend somewhat on 
the three-dimensional structure of the gas and the stellar 
system. The brightnesses of the clusters do not necessarily 
increase monotonically with time due to non-steady gas 
motions, particularly accretion flows in the control simula¬ 
tions. In the feedback simulations, the timescales on which 
the clusters become optically revealed can be substantially 
shorter than the timescales on which the host clouds are 
becoming unbound, and can also vary substantially with 
viewing angle due to anisotropy in the intracluster gas. 
Magnitude-limited surveys may therefore miss substantial 
numbers of clusters simply because of their orientation, and 
it may not be advisable to try to estimate the dispersion 
timescale by attempting to infer the revelation timescale 
observationally. 

(vi) Despite large quantities of damage being done 
to some clouds in terms of gas unbound or expelled, the 
mass or number fractions of stars unbound by feedback are 
very modest and do not correlate very strongly with the 
amounts of gas lost. The cause of this discrepancy is likely 
to be the strong variations in local star formation efficiency, 
which approaches extremely high values in many of the 
clusters or sub clusters. The stars are often situated in 
stellar-dominated potentials in or close to virial equilibrium 
before the onset of feedback, so that the stellar systems are 
largely immune to the destruction of the surrounding clouds. 
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